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Motivated by recent experiments on radiative recombination of two-dimensional electrons in acceptor doped 
GaAs-AlGaAs heterojunctions as well as the success of a harmonic solid model in describing tunneling between 
two-dimensional electron systems, we calculate within the harmonic approximation and the time dependent per- 
turbation theory the line shape of the photoluminescence spectrum corresponding to the recombination of an 
electron with a hole bound to an acceptor atom. The recombination process is modeled as a sudden perturbation 
of the Hamiltonian for the in-plane degrees of freedom of the electron. We include in the perturbation, in addition 
to changes in the equilibrium positions of electrons, changes in the curvatures of the harmonically approximated 
potential. The computed spectra have line shapes similar to that seen in a recent experiment. The spectral width, 
however, is roughly a factor of 3 smaller than that seen in experiment if one assumes a perfect Wigner crystal 
for the initial state state of the system, whereas a simple random disorder model yields a width a factor of 3 too 
large. We speculate on the possible mechanisms that may lead to better quantitative agreement with experiment. 

PACS numbers; 73.20Dx ; 71.45.-d ; 78.20.Ls ; 73.40. Hm 



I. INTRODUCTION 

The properties of non-relativistic electrons in a jellium 
background have been a subject of experimental and the- 
oretical interest. In its quantum mechanical description, 
the Planck's constant h and the electron charge e and 
mass 771, conveniently provide the natural scales for mea- 
suring the largeness/smallness of all other input param- 
eters (such as the average density of electrons) that may 
enter the Hamiltonian of the system. Using the varia- 
tional principle, Wigneiu pointed out that at low elec- 
tron densities a (three dimensional) system of electrons 
with a wavefunction corresponding to a "solid crystalline 
state" would have lower energy than that calculated per- 
turbatively beginning with the non-interacting electron 
"liquid" ground state. Experimental evidence for this 
crystalline "phase" was later presented_for electrons on 
the (two dimensional) surface of Helium.H Variational cal- 
culations by Tanatar and Ceperleya on a two dimensional 
electron system (2DES) supported these observations by 
indicating a phase transition when the electron density 
n was small: (7m)~2 = 37± 5 ( —Vg in units of ;^). In 
recent years, the 2DES has been realized at the interface 
of semiconductor heterostructures. In these experiments, 
the application of a perpendicular magnetic field B pro- 
vides a second parameter, the filling factor v (= ^o^), 
the variation of which has resulted in the experimental 
observation of a rich phase diagrain including integercl 
and fractional quantum Hall statesEl at particular values 
of v. Using the Lindemann melting criterion, the solid- 
liquid phase boundary (iu, the Ts-j^ plane) for the 2DES 
has also been computedE This calculation shows that 
the density at which the liquid to solid transition occurs 
increases monotonically with the magnetic field. It also 
shows that even in the extreme quantum limit (r^ — s- 0) 



the solid phase exists below v = 0.2 . 

Several recent experiments have reported on the ra- 
diative recombination of electrons from the 2DES at the 
junction of a GaAs-AlGaAs single heterostructure with 
acceptor atoms |-a|-Gertain distance from the interface in 
the GaAs layer .Qil3. For certain low values of i^ the pho- 
ton spectrum in these experiments has been interpreted 
as corresponding to the recombination of an electron (of 
the 2DES) in the Wigner solid phase with localized holes 
bound to the acceptor atoms. Accepting this interpreta- 
tion we calculate here the line shape of the spectrum from 
the Wigner solid and compare it with the experimental 
curve.a Our goal is to explicitly investigate whether the 
theoretical consequences of this interpretation agree with 
the actual experimental observations of Ref. ||. 

The model we adopt is as follows. We consider an elec- 
tron solid in which an unoccupied localized state out of 
the plane of the 2DES is available for an electron to tun- 
nel into. One of the electrons - presumably the one closest 
to it - is assumed to tunnel into the available state, leav- 
ing behind a set of electrons that is no longer in an equi- 
librium configuration. This configuration may be repre- 
sented as a distorted solid around some new equilibrium 
configuration, and these distortions may be written as a 
linear combination of the phonon states of the Wigner 
crystal. The squared amplitude for each final state then 
represents the probability of finding the system in some 
definite state of the final Hamiltonian; the difference in 
energy between this state and the initial state is released 
as a photon. Thus one expects a broad photon spectrum 
due to the many possible final phonon states. 

We consider both the limit of a perfect electron solid 
for the initial state, and one in which several charged 
impurities may be located in the acceptor plane. The 
motivation for the latter model is that the experiments 
with which we wish to compare - time resolved photolu- 
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minescence (PL) spectroscopy - measures the spectrum 
only after a waiting time has passed. This not only al- 
lows the electrons to settle down into a quasiequilibrium 
state, but also allows the acceptors that are closest to po- 
sitions of the electrons in the 2DES to absorb electrons 
from it. Thus, after a sufficient waiting perood, the ini- 
tially neutral acceptor atoms, that are meant to probe 
the system, become a source of disorder for it. The dis- 
order associated with the acceptor positions, therefore, 
may affect the PL spectrum. 

Because the holes closest to electrons are eliminated 
quickly, we assume that the holes in our problem are rel- 
atively far from any electrons. In the case of a perfect 
crystalline initial state, we assume that hole,to be located 
near an interstitial position of the lattice.El For a disor- 
dered initial state, we search for locations in the plane 
that are far from any electron and presume these are good 
candidates for locations of holes that may contribute to 
PL long after the initial excitation. In either case, we 
find PL lineshapes that are in good qualitative agree- 
ment with experiment.EI However, for a perfect Wigner 
crystal, we find a PL spectral width approximately a fac- 
tor of 3 narrower than the experimental results, whereas 
our random disorder model yields a lineshape a factor of 
3 broader than experiment. This quantitative disagree- 
ment may arise from correlations in the acceptor posi- 
tions which our random disorder model neglects. 

Our approach is unusual in that we allow not just for 
shifts in the electronic positions, but also for changes in 
the phonon frequencies after the recombination event has 
occurred. This can be a potentially important improve- 
ment, because the removal of an electron effectively re- 
moves a degree of freedom, and may introduce localized 
phonon states in the vicinity of the (charged) acceptor 
atom. Previous approaches to sudden switching prob- 
lems such as this have assumed the phonon spectrum to 
be the same around the init ia l a nd final equilibrium con- 
figurations of the electrons.liilil3 

This paper is organized as follows. In section II we de- 
scribe our model in details and explain how the PL spec- 
trum may be derived from it. Section III is a discussion 
of our results, and we conclude with a summary in sec- 
tion IV. Because we employ novel theoretical techniques 
not readily available in the literature, our appendices are 
fairly complete, which, however, can be skipped unless 
one is interested in following the technical details of our 
calculation. Detail of the calculation are given in these 
Appendices. 



II. THEORETICAL APPROACH 

The harmonic solid model has been successfully used 
previously to study tunneling between two 2DESs in a 
strong perpendicular magnetic field.Eil The quantitative 
success of this nmdel in explaining the experimental I- 
V characteristicsEj has motivated us to adopt a similar 



model for computing the PL spectrum. We therefore 
determine the Hamiltonian /or the x-y degrees of freedom 
of the 2DES (assumed to be in the x-y plane) within 
the harmonic approximation. This approach could be 
considered to be complimentary to the previous Haxtree- 
Fock theory for the phonon shakeup effects on PL.EJ 

The system studied consists of a finite number of elec- 
trons in the x-y plane within a boundary of the shape of 
a parallelogram that is commensurate with a triangular 
lattice. In the model we study there are always some 
electrons that are pinned - they may be in the x-y plane 
and/or below the plane corresponding to the charged ac- 
ceptor atoms. The electron-electron interactions are how- 
ever not purely Coulomb like - this is primarily due to the 
application of periodic boundary conditions which simu- 
lates an infinite system by repeating the finite system at 
integral multiples of its lattice constants. (These lattice 
constants, ai and 02, correspond to the parallelogram 
that serves as a unit cell that contains the entire finite 
system.) In some of the cases studied here the Coulomb 
interaction has been "softened" to the form (r^ -I- z^)~2 
to account for the finite extent of the 2DES in the z di- 
rection. 

The process of electron capture from the 2DES by the 
acceptor atom is modeled to be a sudden perturbation 
of the Hamiltonian corresponding to the x-y degrees of 
freedom of the 2DES. The recombination process is rep- 
resented through the introduction of a fictitious parabolic 
external potential in the final Hamiltonian that confines 
the x-y coordinates of one of the electrons to a position 
To corresponding to the x-y coordinates of the acceptor 
atom that is capturing the electron. The change in the 
z coordinate of this recombining electron is also reflected 
as a sudden perturbation of the corresponding electron- 
electron interactions. The initial {Tii) and final (Tif) 
Hamiltonians, (corresponding to the x-y degrees of free- 
dom) in the presence of a uniform magnetic field in the z 
direction can therefore be written (using the symmetric 
gauge for the vector potential A) as: 

fc=i 

M M 



k=l 1=1 



with 



A{f) = -Biz X r and. 



V{r,z) 




where 



REVT^ 3.0 



Das Sarma Group Preprint, 1997 



R — rtixax + m2a,2 , rrii^ € integers and 
A^\ai X a2\ . 

In the above equations, m* is the effective mass of the 
electron, eo is the dielectric constant of GaAs. fi=c cor- 
responds to the electron (which is necessarily one that is 
not pinned) that is recombining with a hole at ro- The 
limiting prescription in the expression for V{f, z) is for 
removing the singularity in the lattice sum (the summa- 
tion over {R}) which is due, to the long range nature 
of the Coulomb interaction.llj The 6 function in this ex- 
pression is equal to one if the condition in its argument 
is true and zero otherwise - it allows one to represent a 
particle to be interacting with its images in other unit 
cells, but not with itself, z^/^^ is the out of plane separa- 
tion between particles k and / in the initial (final) state. 
A^i(/) is the number of electrons that are dynamical de- 
grees of freedom (in the sense that their coordinates can 
"evolve" ) and M is the total number of electrons includ- 



ing those that are pinned, 
to {4n' ^,^ = 1 to M}, 



Hence, in TL 



*(/): 



in addition 



{ffe, k = N,^f) + 1 to M} 
also serve only as parameters and not dynamical degrees 
of freedom. The choice of these parameters varies with 
the different cases studied here and will be described in 
detail in the next section. Ai = and A/ is chosen to 
be a much larger than the "natural" scale (of the same 
dimension) in this problem: e^n^-^/eo. This choice of 
A/ models the experimental situation in which the elec- 
tron recombines with a screened localized core hole. In 
principle one would like to remove this recombined elec- 
tron as a degree of freedom in the final Hamiltonian, but 
from a calculational point of view it is easier to keep it 
and introduce the last term in Eq. 1, effectively freez- 
ing its motion. The lattice sums in the potential eneiaji 
terms can be evaluated using standard techniquesQO'tll 
(see Appendix A). 

The harmonic approximation to 'Hiu) (denoted by 
TC'^ff.) is developed by expanding the corresponding total 
potential energy V (=the sum of the last two terms in 
Eq. 1) about the (classical) equilibrium configurations. 
The equilibrium configuration is reached by changing the 
coordinates of all the unpinned electrons such that the 
forces on them become zero. The coordinates of electrons 
chosen at the beginning of this "evolution" vary with the 
different cases studied here and will be described in de- 
tail in the next section. The algorithm for this evolution 
is due to Schweigert and Peelers Jlj This algorithm up- 
dates the coordinates by finding the minimum of V when 
expressed to second order in the coordinate increments 
i.e. if V is expressed as a function of the coordinates of 
all the unpinned electrons {qi,l < i < 2N, N = Ni or 
Nf) , then it may be written approximately {V""), to sec- 
ond order in the coordinate increments from the current 
position (g|) as: 



y''(q) = y(q*)-(f*)^(q-q*) + i(q- 



q*)^D*(q-q*), 



with 

ct,T 



(tr 



W(q) 



9q 



and D* = 



dqdq^ 



q=qt 



where q is a column vector whose components are the 
coordinates qi and therefore f becomes a column vec- 
tor whose components are the forces and D becomes 
the (symmetric) dynamical matrix. The lattice sums 
appearing in f A]p4rJ? are again evaluated using stan- 
dard techniquea30llZl (see Appendix A). The updated 
coordinate (q*^^) is determined by minimizing V°'{q). 
However, far away from the equilibrium configuration the 
matrix D* may not be positive definite and hence V^ (q) 
may not have a minimum. Therefore the matrix D* is 
(arbitrarily) changed to D' = D* -I- ry'I where 77* is a 
convergence parameter suitably chosen such that D' is 
positive definite. We have chosen it to be proportional to 
the magnitude of the largest component in f*. In most 
cases this modification of D resulted in a stable evolution. 
However in some cases it failed i.e. D' was nearly singu- 
lar leading to numerical instabilities. In such cases, D' 
was more carefully constructed after diagonalizing D*: if 
D* = MAM"^ (where A is the diagonal matrix consist- 
ing of the eigenvalues of D*) then D' — M/(A)M'^ such 
that /(A) is a diagonal positive definite matrix with the 
function / chosen to satisfy /(A) ^ A as convergence is 
reached. Hence the algorithm for updating the coordi- 
nates becomes: 



,t+i 



{^r 



The above algorithm is iterated until convergence is 
reached - the root mean squared force per electron be- 
coming smaller than a chosen value. This yields the equi- 
librium configuration q^''. 

The process of determining the equilibrium configu- 
ration provides all the necessary parameters (from the 
potential energy terms) for the harmonic approximation 
to the Hamiltonian. In this approximation the potential 
V is replaced by its approximate form V^''(q), with the 
expansion of V°'{q) around the equilibrium configuration 
q'^''. Hence, since f'^'^ = 0, the parameters that enter the 
approximate Hamiltonian are contained in the dynami- 
cal matrix D®'*. This matrix is positive definite since the 
equilibrium configuration reached is stable. There are no 
modes corresponding to neutral equilibrium since global 
translational invariance is broken by the pinned electrons. 
There is no global rotational invariance even in the ab- 
sence of pinning due to the chosen periodic boundary 
conditions. Further, due to the linear variation of the 
vector potential A{r) with r, the approximated initial 
{'Hi) and final Hamiltonians (Ti't) can be written (with 



«(/) 



= ^(q^f/))) as: 



njh 



1 



2m 
1 



(n Bq)^ (n Bq) 

T 



-2*^ 



eq 






(^ -<(/)) 



*(/) 



(2) 
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where 11 is a column vector whose components 11^ are 
the (canonical) momenta corresponding respectively to 
the coordinates % (1 < fc < 27V/) and B is a (anti- 
symmetric) matrix corresponding to the vector poten- 
tial. Note that the harmonically approximated Hamil- 
tonians have an equal number of dynamical degrees of 
freedom {2Nf). For the disordered systems investigated, 
in some cases we chose to randomly pin a small number of 
electrons after the initial equilibrium configuration had 
been found. This is meant to roughiy model the effect of 
pinned in-plane charged impuritiesEil that are not due to 
the charged acceptors. In such cases Nf < Ni. 

Gauge freedom may now be exploited to change A{rk) 
to the form A{fk) = ^Bi^ x {fk - {rk)'i'^) where (ffe)^'' 
corresponds to the equilibrium configuration q^"* of the 
initial Hamiltonian. Further defining #i — q — q^', 
$y = q — q**', and the corresponding canonical momenta 



n, = n, H/ 

written as: 



n B 



q^' - q^ ] , TC^ and TCj may be 






i{f) 



^lu)^if) 



2m* 



m*'n'"' 



B^B -B^ 
B I 



n 



(/) 



with 





= 




+ 


I 

B 


I 
B 




[ qf 1 



(3) 



The harmonic Hamiltonians thus developed can be 
considered to be functions of classical variables or the 
corresponding quantum mechanical operators (with elec- 
trons being spinlcss). In either case, by a linear canonical 
transformation of the phase space variables , they can be 
written as a sum of uncoupled (normal) modes (see Ap- 
pendix B). Hence 






"^U) 



1 

2 






n 



if) 



a 



«(/) 



* 



«(/) 
«(/) 



with 



* 






[c 



«(/)J 






(4) 



(5) 



Since the transformation C^ (C/) is canonical, the com- 
ponents of Hj (H/) are the conjugate momenta corre- 
sponding to the components of the ^i ('4'/) with il; 
(rjy) being a diagonal matrix consisting of the normal 
mode frequencies of Ti'* (Tii). 

The Planck's constant h may now be explicitly intro- 
duced into the Hamiltonians Hi(/) by assuming that 3j(/) 
and ^i{f) are quantum mechanical operators. Hence, 
defining the column vectors consisting of lowering {s^K^f) ) 
and raising {a ,r^) operators: 



Kf)' 



^iif) 



iw) 



I a 
I -a 



* 



(/) 



(6) 



the Hamiltonians may be written as 






(7) 



n 



^Kf)^t(f) 



n 



(/) 



^iU) 



with (using Eqs. 3,5 and 6) 



a/ 

„t 






where T is a matrix of the form 



U V 

V* u* 



I il 
I a 



C/C: 



I I 

ii a 



and w is a vector of the form 



(8) 



, (9) 



(10) 



We now proceed to describe the theory of PL. The ini- 
tial states li/^i) are assumed to be eigenstates oiH^ (with 
^i'lV'i) — Eil'ipi)) whose distribution is determined by the 
temperature P = {ki,T)~^. After the sudden perturba- 
tion, in which the Hamiltonian changes from TCI to Tip 
the initial state is assumed to collapse into a final state 
that is an eigenstate \^f) of H^ (with H'f\ipf) = Ef\^f)) 
with a probability |('0/|i/'i)P- In this process a photon of 
frequency uj is emitted. Assuming conservation of energy 
across the sudden perturbation, hio — Ei — Ef. Hence 
the probability density of the photon frequency V(uj) is 
given by: 



'A 










A* 










(ly 


I 


a 


Cf 


I 


I 




[q-n 


\2n) 


I 


-il 


B 


B 




[^?\ 



V{io) = 



E< 



-/3S. 



(11) 



^e-^^-^|(V'/|V.)P'5(?i-ni?.-i^/)-^) ■ 

i f 

'P{ijj) may be calculated from its Fourier transform V{t) 

1 /■+oc~ 



V{uj) 



2-K 



dter'^^^Vit) , 



(12) 



with 'P{t) given by (inverting Eq. 12 and using Eq. 11) 



V{t) = 



Y^e-^-^ 



(13) 



Rewriting Eq. 13 using the Hamiltonians Ti^ 



H/)- 
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V{t)^ 



E(^^ 



-m, 



\A 



(14) 



ih-^n'}t in-^n'i{t + inf3). 



/ 

Since \ipi) and \ipf) constitute an orthonormal bases, Eq. 
14 may be written in terms of trace (Tr) operations that 
are basis independent: 



Tr 



m^ 



\~in-^n)t^in-^n'i{t + inp) 



Tr 



-pn. 



(15) 



Viiv) is then computed using Eq. 12 . In order to com- 
pare V{oj) to the experimental hne shape of the intensity 
I{(jj) we make the final assumption that: 



I{uj) oc r{uj) 



(16) 



The above relationship (after substituting for 'P{u!) us- 
ing Eq. 11) may be derived using first order time de- 
pendent perturbation theory with a time independent 
Hamiltonian that includes the photon degrees of freedom 
together with the z degree of freedom for the recombin- 
ing electron (see Appendix C). Finally, we note that our 
model for the final state of the recombined electron for- 
mally retains it as a degree of freedom, albeit with a very 
large spring constant A/. In practice this means that the 
spectrum ^(a;) consists of a series of highly separated 
peaks, each corresponding to a different state of the re- 
combined electron plus phonons in the remaining lattice. 
Physically, only the state in which the electron is most 
strongly localized at the location of the core hole - i.e., 
the lowest energy state of the harmonic potential due to 
the Xf in Eq. 1 - is truly relevant. Because of the sep- 
aration of energy scales between the "A modes" and the 
lattice phonon modes, one may easily identify the highest 
energy peak in 'P(w) as the experimentally relevant spec- 
trum. From a computational point of view, one would 
like to eliminate the other peaks, as they are unphysical 
and can consume much cpu time in their computation. 
This can be accomplished by introducing an imaginary 
component i.e., a broadening to the normal mode fre- 
quencies. The precise way in which this is done, along 
with several other practical issues in our computation 
and approximation scheme, is discussed in Appendix D. 
Details of the computation of the traces appearing in Eq. 
15 are given in Appendix E. 



III. RESULTS AND DISCUSSION 

We now present out theoretical results corresponding 
to the experiment of Kukushkin et afl who identify a par- 
ticular "late time" spectrum in their time resolved PL 



spectra as corresponding to the recombination of elec- 
trons from a Wigner solid with holes bound to acceptor 
atoms a certain distance away from the 2DES. For com- 
paring our calculations with their observations, we set 
the electron density n to 5.3 x lO^^/cm^, the filling fac- 
tor ly to 0.1337 (a magnetic field B w 16.4 Tesla), the 
temperature T to 45 mK and the distance between the 
2DES and the acceptor atoms (^o) to 300 A. The effec- 
tive mass parameter m* is set to 0.068?Tie and the dielec- 
tric constant eo is set to 12.8 which correspond to the 
AlGaAS-GaAs heterostructure. We adopt natural units 
defined by n, e^/eo and m* . Therefore the natural length 
scale is {— 1/y/n) « 434 A, the natural frequency scale 



(=V^ 



g2^1.5 



/eo"^* 



1.88 THz and the natural energy 
scale is (= e'^nP'^ /eg) « 2.59 meV. In these units, the 
cyclotron frequency ~ 22.5, h « 0.479 and the inverse 
temperature {(3) ss 668. 

It must b&jiotcd that the experimental late time spec- 
tral "width" EI (the range of photon frequencies for which 
the intensity can be distinguished from being « 0) is 
r-^ 3.0 in the natural units chosen. We study four cases 
here in an attempt to reproduce this width and thereby 
identify the experimentally relevant one. In one of these 
cases we also attempt to give an alternative interpre- 
tation to the experimental observatioia of a "double 
peak" during continuous illumination. As the experi- 
mental spectra have thcjutensity plotted against increas- 
ing photon wavelength,El we show the calculated spectra 
with decreasing photon frequency. The comparison of 
line shape can then be directly made since in the ex- 
periment the spread in photon frequency is a very small 
fraction (~ 2 x 10"'^) of the average frequency. As men- 
tioned previously, the binding energy of the electron to 
the hole, which determines the position of the spectrum 
along the frequency axis, cannot be determined in our 
model and has been set to zero when presenting the re- 
sults. Therefore the calculated peak always appears in 
the negative frequency domain. 



A. Perfect Wigner lattice with Coulomb interactions 

In this case the initial (classical) equilibrium configura- 
tion corresponds to a perfect triangular lattice. Because 
of the translational invariance of a perfect, unpinned 
Wigner crystal, one cannot apply the harmonic approx- 
imation consistently to this case: the sudden change in 
potential will introduce a large coherent motion of the 
center of mass in the subsequent dynamics of the sys- 
tem. This behavior is unphysical, because in practice a 
Wigner crystal will always be pinned by disorder. We 
therefore assume that the electrons at the boundary of 
our supercell are actually pinned, and are not dynamical 
degrees of freedom. This allows the harmonic approxi- 
mation to be applied consistently. We will show below 
that the PL spectrum converges rapidly as the system 
size increases, so that the pinning at the boundary does 
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not affect our final results. 




1 2 



4 5 6 





FIG. 



1 2 3 4 5 6 

1. Equilibrium x-y configuration of the perfect 



Wigner lattice (a) corresponding to the initial Hamiltonian 
Hi (Eq. 1) and the same of the perturbed lattice (b) corre- 
sponding to the final Hamiltonian 7i/. Pinned and unpinned 
electrons are shown by o and + respectively. The figures cor- 
respond to the case when the number of unpinned electrons 
Ni = Nf = 9 with the total number of electrons M = 16. 
The "central" (unpinned) electron is shown by *. The axes 
are labeled in units of the natural length scale of the problem; 
(total electron density)" 



-1/2 



Here, all the electrons interact through the "ideal" 
Coulomb interaction. The initial configuration can there- 
fore be directly constructed by laying down the electrons 
at lattice positions or arrived at by evolving an initial 
configuration in which only the pinned electrons are laid 
down at lattice positions with the rest of the electrons be- 
ing randomly placed. The final configuration, obtained 
through evolution beginning with the perfect Wigner 
lattice, corresponds to confining the "central" electron 
(see Fig. 1) to one particular point on the boundary of 
its Wigner-Seitz cell (corresponding to the prefect lat^ 
tice). This choice is motivated by the interpretationEl 
that the late time PL spectrum corresponds to recom- 
bination events in which the distance between the hole 
(bound to the acceptor) and the recombining electron is 
maximal. It must be noted that to determine the final 
equilibrium configuration it is not appropriate to begin 
with a random configuration of unpinned electrons since 
then the equilibrium configuration reached may involve 
exchanges of electrons with respect to the perfect lattice. 



In this work, we assume that the initial state [Fig. 1(a)] 
may be regarded as a deformation of a classical equilib- 
rium that is closest in configuration to this state [Fig. 
1(b)], and that the subsequent motion of the electrons is 
due to the vibrations around the latter state. While final 
states in which electrons are exchanged are in principle 
relevant, such exchanges in practice contribute little to 
the PL spectrumllj 
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FIG. 2. Calculated spectra from transitions similar to 
that in Fig. 1 (which corresponds to the number of unpinned 
electrons (Nf) being equal to nine). For each value of Nf, 
fiat regions correspond to the zero of Vamiu). %err is the 
increase in calculated spectral width due to approximations 
made to get a continuous curve for the finite Nf. The spectral 
width saturates (around Nf — 25) to ~ tt of that observed 
by Kukushkin et al. However, a qualitative feature of the 
experimental spectrum is present in the theoretical results: a 
faster rising edge as compared to the falling edge. 

We can therefore fix the parameters in 'Hi(f) (Eq. 1): 
Ni = Nf ^ P^, M = (P+l)^ (Fig. 1 corresponds to the 



'f 



k.l 



case when P = 3), Zj ' =0 (corresponding to all electrons 
initially being in the x-y plane) and z^' — zo((5;,c — ^k.c) 
(corresponding to all the electrons finally being in the 
x-y plane except for the central electron c being below 
the plane at the position of the acceptor atoms) . tq f^nd 
other parameters corresponding to the x-y position of the 
pinned electrons {fi, i — Afi(/) + 1 to M) can be 
inferred from Fig. 1. We have studied the cases with 
P — [3,5,7,9,11,13]. In each of these cases we have 
chosen the parameter Fg (a broadening parameter intro- 
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duced to make the phonon density of states continuous; 
see Appendix D) to be equal to the smaUest phonon fre- 
quency in 7i,'\ These values (in units of ^J e^n^-'^ j ti^m*^ 
are [7.59, 4.49, 3.08, 2.30, 1.82, 1.49] x 10"^ (corresponding 
to the cases with P = [3,5,7,9,11,13] respectively). It 
must be noted that all the electrons including the pinned 
ones are included in the computation of the density n. 

Fig. 2 shows the calculated spectra. The error in the 
spectral width is due to the approximation of replacing 
the delta function in Eq.ll by a Gaussian of width Fq. 
Note that the spectral width seems to have saturated 
beyond Nf fa 2p but it is only ~ 1/3 of that experimen- 
tally observed.El A qualitative feature of the experimental 
spectrum is however reproduced: a faster rising edge as 
compared to the falling edge. 
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B. Perfect Wigner lattice with softened Coulomb 
interactions 

We now soften all the in-plane electron interactions in 
an attempt to see if it results in a broader spectrum than 
the previous case. The softened interaction is of the form 
1/ y/r'^ + Zg. Zq is chosen to be 150 A - a significant frac- 
tion of the inter-electron distance. This form of the in- 
teractions is motivated by the fact that the electrons are 
not strictly confined to two dimensions at the interface 
of the GaAs-AlGaAs heterostructure, there is a finite ex- 
tent in the z direction of the wavefunction confining the 
2DES. This particularly simple form of the Coulomb in- 
teraction incorporating the finite width of the electron 
layer has been quantitatively successful in obtaining the 
fractiopal quantum Hall excitation gap energies in the 
2DES,B2l which motivates us to apply this softened inter- 
action in the Wigner crystal PL calculation. The proce- 
dure to obtain the initial and final equilibrium configu- 
rations remains the same as before. These configurations 
are similar to that shown in Fig. 1 - there being a "wall" 
of pinned electrons on the boundary. Only one system 
size is studied with N,(f) = 7^, M = 8^, z'l'^ = zo(l-4,i) 



and 



k,l 



'«(/) 



zoil - Sk,i) + {zo - zo)\Sk,c- Si^c\ where c again 



corresponds to the central electron. The parameters tq 
and { fi, i — Ni(^ff -I- 1 to M} can be inferred from Fig. 
1. Again, Fq has been chosen to be the smallest phonon 
frequency of H ■' (= 2.52 x lO^^^e^ni-^/eom*). 

Fig. 3 shows the calculated spectrum. The width is 
slightly reduced (still roughly 1/3 of the experimental 
value) as compared to the case with ideal Coulomb inter- 
actions. This can be attributed to a faster falling rcdge 
making it less similar to the experimental spectrum.El We 
conclude that only softening the Coulomb interaction by 
the finite layer thickness of the 2DES cannot account 
for the experimental spectral width. This is understand- 
able, because the spectral width of the PL spectrum is 
primarily determined by the physics of the recombination 
process at the acceptor sites and not by the details of the 
electron-electron interaction. 
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FIG. 3. Calculated spectrum from transition similar to 
that in Fig. 1 (which corresponds to the number of unpinned 
electrons Nf being equal to nine) with softened Coulomb in- 
teraction. For comparison the spectrum corresponding to 
ideal Coulomb electron interaction is included. It has been 
shifted to the left to make the photon energy corresponding 
to the difference in the ground state energies coinside with 
the same point of the first curve. Softening the electron in- 
teractions results in the spectrum being less similar to the 
experimental curve: the width is slightly reduced ( continues 
to be ~ 1/3 of the experimental value.) and it has a faster 
falling edge. 



C. Perfect Wigner lattice Avith recombination 
position averaging 

Since the "softened" Coulomb interaction does not lead 
to any improvement in the agreement between theory and 
experiment, we now revert back to "ideal" Coulomb in- 
teraction for in-plane electrons and consider averaging 
spectra corresponding to different x-y positions of the re- 
combination center (= acceptor atom position) relative 
to the initial position of the recombining electron. The 
averaging is carried out assuming that the x-y distribu- 
tion of acceptor atoms is uniform. This results in a uni- 
form distribution for the position of the recombination 
center within the Wigner-Seitz cell of the recombining 
electron. It models the situation for which the PL is not 
time resolved Jallj so that no particular final position of 
the recombining electron is favored. Symmetries can be 
exploited to reduce the actual region of the Wigner-Seitz 
cell that needs to be explored. For configurations of the 
type shown in Fig. 1 (a), due to the presence of the 
"wall" of pinned electrons on the boundary, there are 
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only two reflection symmetries (each along a diagonal of 
the system) for the recombination positions correspond- 
ing to the central electron. Therefore, to get the average 
spectrum, 1/4 the area of the Wigner-Seitz cell needs to 
be explored. However, if there was a hexagonal boundary 
of pinned electrons around the recombining electron then 
there would be six-fold together with reflection symme- 
tries along the diagonals which would reduce the region 
needed to be explored to 1/12 the area of the Wigner- 
Seitz cell. We therefore consider such a configuration 
(Fig. 4). 



they are included in the computation of the total energy 
- one third of which is the actual energy of the hexagonal 
system. During evolution, the forces and the dynami- 
cal matrix need to be computed corresponding to only 
the unpinned electrons within any one hexagon (inter- 
actions between all inequivalent electrons are taken into 
account). By symmetry they are the same for the corre- 
sponding unpinned electrons in the other two hexagons. 
The displacements are therefore computed for only one 
set of unpinned electrons but it is applied to all the un- 
pinned electrons thus preserving the equivalence 
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FIG. 4. The figure shows the initial equilibrium configuration of the system with a hexagonal wall of pinned electrons. 
Electrons with the same number are considered equivalent. This, in addition to application of periodic boundary conditions 
about the parallelogram, results in the hexagonal system being repeated periodically - three copies of which are embedded in the 
parallelogram. Therefore the Wigner-Seitz cell of the "central" recombining electron (numbered 1) has all the symmetries of a 
triangular lattice. This reduces the region of different recombination positions to be explored to yrj the area of the Wigner-Seitz 
cell - the triangle within the hexagon around the electron numbered 1. 



Fig. 4 shows the construction of the system with a 
hexagonal wall of pinned electrons. Three such copies 
are needed for embedding into a parallelogram - then the 
usual periodic boundary conditions about the parallel- 
ogram results in the periodic repetition of the hexago- 
nal system. The initial equilibrium configuration is di- 
rectly constructed as it corresponds to a perfect triangu- 
lar lattice. The final equilibrium configuration is arrived 
at through evolution beginning with this configuration. 
However, as electrons with the same number (see Fig. 
4) are considered equivalent, the corresponding "self in- 
teraction" terms are neglected during the computation 
of the forces and the dynamical matrix. Nevertheless, 



of the electrons with same number through the evolution. 
The parameters in 7ii(/) used are: Ni = Nf = 61, M = 

75, Zj ' =0, z/ — ZQ(5k,i—5i^i). The parameters r^, i = 
Nii^f-f -|- 1 to M can be inferred from Fig. 4. 

PL spectra are computed for different recombination 
positions (the parameter ro in Eq. 1) of the electron 
numbered 1 (see Fig. 4). Each of the spectra is computed 
with Fq being equal to the smallest phonon frequency of 
Hf: = 3.46xl0-2y/e2ni-5/eoTO*. Fig. 5 shows then vari- 
ation for recombination positions distributed uniformly 
along the edge of the Wigner-Seitz cell. Although the 
overall normalization increases rapidly with decreasing 
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distance of the recombination center from the initial po- 
sition of electron 1 , the position of the peak shifts only 
slightly to the left. Hence the average of these spectra 
continues to have the typical width o£,each curve: ^ 1/3 
of that observed in time resolved PLH 

The average PL spectrum - using 127 distinct points 
uniformly distributed in the Wigner-Seitz cell - is shown 
as an inset to Fig. 5. The spectrum is relatively nar- 
row, and is in fact dominated by PL from recombination 
events in which the x-y motion of the electron is quite 
small. This is expected as it reflects the fact that the 
overlap between the recombining electron's initial and 
final states falls off rapidly as a function of its displace- 
ment. The spectrum is much narrower than seen in ex- 
periment, and does not refle<:jt||-the experimentally ob- 
served double peak structurcHllO The latter has been 
interpreted as evidence of liquid-solid jaaexistence due 
to disorder and/or finite temperatures jS'tHlEJ, and our 
present results implicitly agree with this interpretation 
- a model of a perfect uniform solid as used in our work 
cannot explain the structure seen in continuous wave ex- 
periments. 
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FIG. 5. Spectra for recombination positions uniformly dis- 
tributed along the edge of the Wigner-Seitz cell (see Fig. 4). 
They are shown in the order of decreasing distance from the 
initial position of the recombining electron. The final spec- 
trum is the average of these spectra. It continues to have a 
width that is tt of that in experiment. * The inset shows the av- 
erage spectrum corresponding to points distributed uniformly 
over the entire Wigner-Seitz cell (its axes are in the same units 
as the figure). It is dominated by the spectra corresponding 
to recombination positions close to the initial position of the 
recombining electron. 



D. Wigner lattice with disorder averaging 

We finally consider a disordered system (again with 
ideal Coulomb interactions) in which the disorder is due 
to pinned electrons below the x-y plane corresponding to 
charged acceptor atoms i.e. we assume that the spec- 
trum of interest corresponds to a time when most of the 
recombination events have occurred so that all the ac- 
ceptor atoms are charged in equilibrium. We consider 
a system size of 64 (initially) in-plane electrons. Unlike 
the previous cases, only these electrons are included in 
the computation of the density n. As the experimental 
density of acceptor atomsEl is 5 x 10^/cm^ the number 
of charged acceptor atoms corresponding to 64 in-plane 
electrons must be very nearly 6. We therefore have 5 
pinned electrons below the x-y plane corresponding to 
these atoms - the sixth would correspond to the recom- 
bination event for which the spectrum is computed. 

In the absence any concrete information about the 
correlations between different acceptor atom coordinates 
that exist due to the process of 5 doping we have assumed 
that each acceptor atom's x-y position is uniformly ran- 
domly distributed and is independent of any other accep- 
tor atom's position. Therefore the 5 pinned electrons are 
independently placed randomly in the system. In princi- 
ple, one would like to find the point farthest from all the 
electrons for a given disorder realization, and, since ac- 
ceptors located near such points are least likely to have a 
recombination event, assume that the late-time PL spec- 
trum is dominated by an acceptor at this point. In prac- 
tice, rather than generate a large number of disorder real- 
izations, we choose just one, construct the Wigner-Seitz 
cells for the disordered system, and choose the corners of 
these cells as candidates for late-time PL recombination 
events. We believe this procedure will produce spectra 
qualitatively and even quantitatively close to that found 
by direct disorder-averaging, and is numerically much less 
time consuming. Essentially, we are assuming the system 
self-averages. To model the likelihood that a particular 
corner-configuration is likely to be available in the late- 
time spectroscopy, we additionally introduce a weighting 
factor for each corner. We do this by constructing the 
Voronoi cell around the lattice of corners, and set the 
weight for the corner to be proportional to its dual cell's 
area. Thus, acceptors that are particularly far from elec- 
trons are more likely to be available for PL events after 
long-time, and are given somewhat larger weights. 

Two cases are studied here: Case (a) - None of the 64 
in-plane electrons are pinned. Case (b) - Two of the in- 
plane electrons that are maximally distant (as measured 
in the x-y plane) from the pinned electrons are addition- 
ally pinned after obtaining the initial equilibrium config- 
uration (electrons numbered 13 and 47 in Fig. 6). The 
latter case,is chosen to study what effects in-plane pinned 
electronstd could have on the PL. As will be seen below, 
their effect is quite small. The parameters of 7ii(/) (Eq. 
1) can therefore be set to: Ni — 64, Nf — 64 for case (a) 
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and 62 for case (b), Af = 69, 

k,l 



k,l 



Zo{Sk-. 



3i,>65- 



where c corresponds 



to the recombining electron. The parameters fi, i — 
^iU) + 1 to Af can be inferred from Fig. 6. 

The initial equilibrium configuration is reached 
through evolution beginning with a random distribution 
of the in-plane electrons. Beginning with this configu- 
ration, the final equilibrium configuration is reached in 
ten steps - each step localizing the recombining electron 
(through the choice of the parameter tq in Eq. 1) closer to 
its final position which is chosen to be one of the corners 
of the Voronoi cell that initially contained the electron 
(see Fig. 6). This procedure allows us to find the closest 
minimum energy state in configuration space. 



der the corresponding spectrum {— Varn{t = 0)). Fig. 
7 shows the normalized distributions for the two cases 
studied here. It may be seen that they have a width 
of ~ 1 ^ e^v}-^ I e^m* . The insets show the correspond- 
ing distributions for a second disorder realization - they 
are similar in shape but have a somewhat smaller width: 
~ 5y/e'^n^-^/eom*. We therefore expect that the average 
spectra corresponding to the first disorder realization to 
be roughly similar to the true disorder averaged spectra. 
The individual spectra in the average spectrum corre- 
sponding to the two cases studied here have been com- 
puted using Fq = 3.80 X lO'^^/e^n^-^/eom* - the smallest 
non-zero phonon frequency of the perfect unpinned lat- 
tice with 64 electrons. Fig. 8 shows these average 








5 







FIG. 6. The figure shows the initial equilibrium configuration of unpinned electrons (marked by -I-) and a particular final 
configuration of electrons (marked by *) corresponding to the recombination of the electron numbered 1 with a bound hole at a 
corner of the Voronoi cell that contained this electron. For easy visualization of the Voronoi cells, the initial configuration has 
been shown to be periodically repeated. Electrons that have already recombined are represented as pinned electrons (marked 
by o and numbered 65 through 69) that are independently and randomly placed in the system. They are not considered in the 
Voronoi construction. 



It appears that the width of the average spectra cor- 
responding to the two cases studied is mainly due to the 
wide distribution [Diajp]) of the position of the peak 
{ijjp) in the individual spectra that constitute the av- 
erage. This distribution is constructed by binning the 
peak positions into intervals in frequency. All peaks 
within a bin do not contribute equally to the distribution 
height - each contributes in proportion to the product of 
the corresponding Voronoi corner area and the area un- 



spectra. 



As expected they are have a very large (> 
7 ^ e^n?--^ / eom*) width: - QJe'^n^-^ /torn* . This is ~ 3 
times the experimental valueB The extra broadening seen 
in Fig. 8 over that of Fig. 7 is due to the shakeup of 
phonons. 

We speculate that the disagreement in these results 
and those found in experiment may be due to an over- 
estimate of the disorder strength assumed in our uncor- 
related random disorder model. In particular, the as- 
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sumption that nearly all the holes on acceptor atoms 
have recombined with electrons would be true only af- 
ter very long waiting times, and it is undLear that this 
limit is actuaily achieved in experiment .0 Subsequent 
measurementaD have suggested that it may indeed be 
the case that longer waiting times are needed. Clearly, 
one could "tune" the level of disorder in our model to 
obtain the experimentally observed linewidth. However, 
this would require a painstaking and somewhat artificial 
"fine-tuning" of our model. An independent measure- 
ment of the density of charged acceptors in these systems 
in the late time PL would greatly facilitate a quantitative 
comparison of this model with experiment. In this con- 
text it is worthwhile to point out that in general modula- 
tion 5-doped heterostructures are knownE3 to have sub- 
stantial correlation among the impurity sites, and an un- 
correlated random disorder model overestimates the dis- 
order strength producing lower electron mobilities than 
experimentally observed. Any such correlation among 
the acceptor sites in our problem would reduce the spec- 
tral width of our calculated PL spectra, bringing experi- 
ment and theory into closer quantitative agreement. 
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FIG. 7. The figure shows the distribution of peak posi- 
tions (extended flat regions correspond to its zero value) in 
the spectra corresponding to the disorder reahzation shown in 
Fig. 6 for the two cases: (a) with no in-plane pinning centers 
and (b) with two in-plane pinning centers. The broad distri- 
butions imply that the corresponding average spectra must 
be at least as wide. The insets (whose axes are labeled in the 
same units as the figure) show the corresponding distributions 
for a second disorder realization. They are similar those of the 
first but have a slightly smaller width. This suggests that the 
two average spectra corresponding to the first disorder real- 
ization must be roughly similar to the true disorder averaged 
spectra. 
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FIG. 8. The figure shows the average spectra correspond- 
ing to the disorder realization shown in Fig. 6. Case (a) corre- 
sponds to all in-plane electrons being unpinned whereas case 
(b) corresponds to there being two additional in-plane pinning 
centers. As expected (see Fig. 7) the average spectra are very 
wide: roughly three times that observed by Kukushkin et al.^ 



IV. CONCLUSION 

The line shape of the PL spectra due to electron re- 
combination from a finite and pinned two dimensional 
Wigner crystal to a hole bound to an acceptor atom, 
computed using the harmonic approximation and first 
order time dependent perturbation theory, is similar to 
that seen in experiment -it has a faster rising as compared 
to the falling edge. However, the width of the calculated 
spectra, for recombination events beginning with the per- 
fect lattice configuration of electrons, is only the ~ 1/3 of 
that observed in experiment E With the initial configura- 
tion of electrons being disordered due to charged accep- 
tor atoms corresponding to already recombined electrons, 
the spectral width is about three times the experimental 
value. We have speculated that considering lesser disor- 
dered configurations and recombination events different 
from those studied here may result in better agreement 
with experiment. Since disorder in our model is due only 
to charged acceptor atoms we speculate that any corre- 
lations in their positions (that have been assumed to be 
zero here) may also reduce the corresponding computed 
spectral width. Because our perfect crystal and random 
disorder calculations give results factors of three smaller 
and larger than the experimental PL width respectively, 
one could perhaps get quantitative agreement with the 
experimental spectra (we emphasize that our theoretical 
results are in good qualitative agreement with the exper- 
imental PL spectra) by using an adjustable correlated 
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disorder model, but we feel that, in the absence of any 
concrete information about the nature of disorder, this is 
not a particularly meaningful exercise. 
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APPENDIX A: EWALDS SUMS 

The lattice sums appearing in the potential energy 
terms of Eq. 1 can be obtained as a special case of the 
following sum for E: 



E{r,q,z,p,d) = ^ 



Jq- R 



R 



e{f+ze^ = 0) 
\r + ze, + R\P \r + ze,\P 



0{p<d,q^O) f d'^r'e^'i-'' 



A 



\r + zcz + r 



TTTv^ (Al) 



with i? = > micii, rrii G integers, di ■ Bz = OV i ■ 
i = 1 

In the above formula a^ are the primitive lattice vectors 
of the d dimensional direct space lattice, and A is the unit 
cell volume, r is assumed to lie within the region spanned 
by the unit cells touching the origin of the direct lattice. 
Similarly, q is assumed to lie within the region spanned 
by the unit cells touching the origin of the reciprocal 
lattice, the unit cell of which is defined by the vectors bi 
that are determined by the following: 

di ■ bj — 2Tr6ij . 

The unit vector Bz in Eq. Al does not belong to the 
space spanned by di and lies along an additional dimen- 
sion "perpendicular" to the direct space lattice, p may 
be any positive real number. The function is equal to 
one if all the conditions in its argument are true and zero 
otherwise. In the Madelung sum, we are interested in 
the interaction energy of a particle with all other parti- 
cles, including its images in other unit cells. This is most 
conveniently handled by a sum over the Bravais lattice 
vectors (first term in Eq. Al), but then one must remove 
the interaction of a particle with itself (second term in 
Eq. Al). For long ranee interactions {p > d), such sums 
diverge unless an interaction with a neutralizing back- 
ground is introduced. This is the meaning of the last 
term in Eq. Al. V that appears in Eq. 1 may now be 
written as: 



U(f,z) = -i?(f,0,z,l,2). 

Each term of the lattice sum (LS) (the summation in 
Eq. Al) for E can be written as a sum of a rapidly de- 
creasing function /i (R) and another function /2 (R) that 
is slowly varying. This division is achieved by multiply- 
ing the terms by the following expression for unity in 
terms of gamma functions: 



1 = 



Tp 
1 



rfP 2 2\ I /P 2 2\ 



where 



If+ze^ + i?! , r(a,6) 



e-H^-^dt 



ra = r(a,o) 



7(a,5) = ra-r(a,6) , 



and e is a suitably chosen (inverse length scale) conver- 
gence parameter the choice of which is explained later 
in this section. Hence the functions /i and J2 may be 
identified as: 



MR) 



Tp'x-Pe^'l-^'ril^e'x') 



(A2) 



with the LS given by 

L5 = ^/i(i?) + ^/2(i?) 
R R 



(A3) 



Due to the functional form of the gamma function 
T{a,b), the terms of the first summation in Eq. A3 die 
off rapidly with increasing e|.R| i.e. for e\R\ ^ 1, 



\hiR)\^TphP{e\R\)-'e-i'\R\) 
1 



(A4) 



Therefore this summation is carried out directly over {i?} 
(beginning with i? = 0) until the estimated relative error 
due to the neglected terms is smaller than a chosen value. 
Due to the functional form of the gamma function 
7(a, 6), the terms of the second summation in Eq. A3 
vary slowly with increasing e|-R|. Therefore their summa- 
tion is better performed using the following identity that 
transforms the summation over {i?} of /2 to a summa- 
tion over reciprocal lattice vectors ({G}) of the Fourier 
transform of /2 (/2): 



R 



R 



f2{r')W)d'r' , Sp{f') =J2Sir' 

R 

1 



R) 



(27r) 



f2{k)5p{-k)d''k , 
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where f2{k) ^ I f2{r")e^'' ' "" d'^r' 



6p{k) = / Spiry'^-'^'d'^r' = i^^<5(fc ^ G) 



A 



It may be rewritten, separating the possibly singular 
terms (corresponding to R = and G = ) and evalu- 
ating them as appropriate limiting cases (of f+ ze^ -^ 
and q ^ respectively) as: 



G 



with G = 2, ^TT-A , rrii G integers, so that 
i = 1 

R G 

f2{k) is computed to be: 



h{k) 
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2e 
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2e 



,eV), (A5) 



with y = |G — (fl 



It must be noted that T(a, 6, c) < T(a, 6, 0) = 
6^°r(a, 6). Therefore, again due to the functional form 
of the gamma function r(a, 6), the terms in the second 

IGI 
summation of Eq. A5 die off rapidly with increasing ■U-'- 

i.e. for J^ > 1 



l/2(-G)| < f (\G\ 



2e 



(A6) 



Therefore the second summation in Eq. A3 is trans- 
formed using A5 and is carried out over {G} (beginning 
with G = 0) until the estimated relative error due to the 
neglected terms is smaller than a chosen value. 

E can now be written (using Eqs. A1,A2,A3 and A5) 
as: 

E{r,q,z,p,d) = 
e 



f ^ {e\r + ze, + R\)P 



\r + ze 



R 



f 



r~^eP 
I y,^(G-q)-r^(d_v(\G_£^ 2^2^ 

A{elV^f^ '^ \^^) 

'^ J \f+ze,+f"\P ■ 



Eif,q,z,p,d) = 



(A7) 



+ - 



TphP 

I p e 
J, 

Aie/^Y 



\f+ze^ + R\ 



J^ {e\r + ze, + R\)P 

R^O 



ye 



-tq ■ r 



^ ^t{G-q)-r^^d-p^ 
G^O 



'2e 



,e^z^) 



with 



r(|,.M.' 



2 
P 



e\r + ze^\P 



y = 



T( d~P ( \~q\ 



-T( 



p — d 



,e'z\Q) 



if r -f ze^ 7^ 
if r + ze^ = 

, e^z^) if p > d or (f 7^ 
if p < d and q = 



Equation AT can be used to evaluate E for all values 
of its arguments. First and second derivatives of E (with 
respect to r) that are needed to compute the forces and 
the dynamical matrix are also obtained by performing 
the requisite operations on the above expression. For 
equally rapid convergence of the two summations in this 
equation, e must be given by: 

(A8) 



A 



1/d 



With this choice of e, the magnitudes of the largest term 
neglected (|(5T|), if TV of the largest terms are included 
in each of the two summations, become equal - and may 
be estimated (using Eqs. A4 and A6) to be: 



\5T\ 



P 
1 



^y^r^-^' 



e 



A 



r+l 



.1/d 



From previous work,EZI T(a,5, c) can be written in 
terms of Gamma functions for a being any odd multi- 
ple of ^ (and b ^ 0) using: 

T(a, b,c)^^ e"(^ + ^) + aT{a - 1, b, c) + cT{a - 2, 6, c 
T(i,6,c) = ^[a + /3] , T(-i,6,c) = ^[a-/?] , 
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with 



a = e 



-2V6c 



'l-[0(6>c)-0(&<c)]7(i 1%/6-V^I) 



(3 = e'^^^ri\,Vb + V-c). 

Finally, an accurate evaluation of Gamma functions {rel- 
ative error p^i-lJ^"^^) is done using the ideas outlined by 
C. Lanczosf* 



APPENDIX B: NORMAL MODES OF A 
GENERAL HARMONIC HAMILTONIAN 

This section shows the procedure for transforming a 
harmonically approximated Hamiltonian of the form 



n'^^e" 



2m* L J 



-B 



B^ 
I 



n 



into a summation over uncoupled (normal) modes by a 
linear canonical transformation of the phase space vari- 
ables [^-^n-'"] (Eqs. 4 and 5). The only requirement 
for this procedure to apply is that D be a symmetric 
positive-definite matrix. This is the case here since D is 
the dynamical matrix corresponding to a stable equilib- 
rium. (All matrices that appear in this section are real.) 
The transformation is carried out in three steps which 
together constitute a canonical transformation. This is 
verified by simultaneously determining the equations of 
motion for the transformed variables and showing that 
they may be derived from the Hamiltonian in the same 
variables. The equations of motion in the current vari- 
ables read as: 



n 



1 

m 



-B I 

m*D-B^B B^ 



n 



The procedure is similar to the standard technique^Ea but 
has to be extended to a multiparticle generalization. 

The first of the three stpes in the transformation is 
in itself a canonical transformation that diagonalizes the 
dynamical matrix D. This is accomplished by an orthog- 
onal matrix N (i.e. N-^ = N~^) the existence of which 
follows from the real symmetric structure of D. Since 
the eigenvalues of D are all positive (as it is a positive 
definite matrix) the diagonalized result may be written 
as: TO*D^ = N-^DN where D^ is a diagonal matrix 
with positive elements of the dimension of frequency. The 
transformed variables may therefore be written as: 



Hi 



(m*D^)^/2N^ 



(Bl) 







(m*D^)-^'^N' 



1/2. 



n 
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in terms of which the Hamiltonian and the equations of 
motion read: 



n'' = e" 



*i Hi 



ill 



D, 



D, 



BjDjiB. 



B^D-iB. 






B^ 



*i 
Hi 



*i 
Hi 



where B^ is the transformed B matrix that also has 
the dimensions of frequency and is given by: m*'Q^ = 
DJ^/^N'^BNO^^ It may be seen that if B = 
(=> Bt^ = 0) then Eq. Bl would suffice to diagonalize 
Ti^ into the required form. 

The next transformation trivially diagonalizes the 
Hamiltonian. However, as the equations of motion in 
the new variables do not follow from the Hamiltonian in 
the same variables, it is not a canonical transformation. 
The transformation is given by: 



*2 

Ha 



D 



1/2 







D 



-1/2 



B^ D 



1/2 



*1 

Hi 



(B2) 



In terms of the new variables the Hamiltonian and the 
equations of motion now read as: 



*2^n2^ 



n 


'' = . 


ni 1 1 
^ 2 


"*2 







.n2. 




.-D^ ] 



M - B„ 





I 






*2" 




I 




n2 




*2 






n2 









Q 



*2 

n2 



where B^ = TiZ^'^B^Y>]J^ . 

The final transformation, which together with the pre- 
vious one results in a canonical transformation, involves 
the "block off-diagonalization" of the matrix Q by an or- 
thogonal matrix M {i.e. M* — M~^) the existence of 
which follows from the real antisymmetric structure of 
Q. It must be noted that Q is an even-dimensional ma- 
trix (27V/ X 2Nf) which may be shown to be non-singular 
since all the diagonal elements of D^^ are non-zero. M is 
constructed using the eigenvectors of another symmetric 
matrix Q: 



Q 



Q 
Q 



It can be shown that the eigenvalues of Q are neces- 
sarily non-zero (due to Q being nonsingular) , degener- 
ate by an even number and that corresponding to ev- 
ery eigenvalue there is another with equal magnitude 
and opposite sign. This can be seen from the construc- 
tion of the eigenvectors: if there is an eigenvector of the 
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form [vi"^, V2^]"^ with eigenvalue wi (vi and V2 are col- 
umn vectors with 2Nf components), then the vectors 
[v2'^, — vi^]"^, [v2'^,vi^]^ and [— vi^,— V2^]^ are also 
eigenvectors with eigenvalues oji, —lui and — o^i. It may 
be shown that the four vectors are mutually orthogonal 
and hence linearly independent. From the mutual or- 
thogonality of all eigenvectors of Q constructed in the 
above fashion (for eigenvalues that are more than two 
fold degenerate, the corresponding eigenvectors may need 
an explicit orthogonalization procedure) it can be shown 
that the vectors Vi*, V2' (i = 1 to Nf) are mutually or- 
thogonal. They are chosen from the set of eigenvectors 
of Q that have coi > and are explicitly normalized. The 
matrix M is then constructed: 



M=[vi\vi2, ,V1^/,V2\V2' 



,V2 



Nf 



It may be verified that M is orthogonal and that it "block 
off-diagonahzes" Q: 



M^QM 



n 
n 



where H is a diagonal matrix with the diagonal elements 
being the (positive) phonon frequencies uji. The final 
transformation is given by: 



* 



0-1/2 



[M-] 



*2 
II2 



(B3) 



In terms of the variables ^ and H the the Hamiltonian 
and the equations of motion read: 



n 


" = 6 


™ + 1 r^^H^- 




no 
n 




* 


r*i 




r ^] 




[*] 




[I]- 




n 




_ 3 _ 











As the equations of motion for ^ and H follow from the 
Hamiltonian in the same variables, the combined trans- 
formation C represented by Eqs. Bl, B2 and B3 must 
be canonical. In other words the transformation C given 
by: 



n 1/2 

f|-l/2 



[M-] 



D 



1/2 



(m*D^)"'/" 



J-'cJ -L'^ ■'-'1 







1/2 



N^ 



APPENDIX C: PHOTOLUMINESCENCE 

FORMULA FROM TIME DEPENDENT 

PERTURBATION THEORY 

In this section we postulate a time independent Hamil- 
tonian Ti" which, on the application of first-order time 
dependent perturbation theory, yields the PL formula 
used (Eqs. 16 and 11). (In this section all bold faced 
vectors have three components (x, y and z) whereas oth- 
ers have only two (x and y): v = {v,Vz).) 

The Hamiltonian 7i° may be assumed to be an ap- 
proximation to another Hamiltonian Ti. that has been 
postulated previouslyES to include the (first quantized) 
degrees of freedom for (spinless) electrons (corresponding 
here to unpinned electrons, Nf in number) with the elec- 
tromagnetic fields considered to be the sum of two fields 
- internal fields that are due to the electronic degrees of 
freedom and external fields (here the uniform magnetic 
field in the z direction and the electric fields due to pinned 
electrons and any required confinement potentials) that 
are independently specified. The internal fields are con- 
sidered dynamical - their energy therefore included in Ti. 
These fields are represented through a four vector po- 
tential {Ap,(J>p) in the Coulomb gauge: V • A = 0. It 
may be seen that in this gauge A^ is dynamical (and 
therefore quantized) whereas <^^ is not - it is completely 
determined by the electronic degrees of freedom as the 
total Coulomb potential. Its contribution to the field en- 
ergy results is the Coulomb interaction between the elec- 
trons. The external fields are considered non-dynamical 
and therefore do not contribute to the field energy. They 
are specified through another four vector potential (A, 0) 
that is a function of the electronic degrees of freedom and 
whose gauge is chosen independently. Both the types of 
fields are minimally coupled to all the electronic degrees 
of freedom. This gives: 



Ti = Tio + Til ; 'Ho = Ti^ + TYq , 



(CI) 



where 



1 

2m 



e -• 1 -^ ■r-^ 

Pfc - -A(rfc) + e(/)(rfc) + > 



l>k 



rfc - r; I 



U / ^ k !(■ ^ k,r ' 

k, r 

-e 



^^ = E 



-AP(rfc) 



Pfc A(ffc) 

c 



2mc 



AP{rk) 



may be verified to be canonical, i.e. it satisfies: 



I 
I 



C therefore transforms the Hamiltonian into uncoupled 
normal modes as required (Eqs. 4 and 5). 



with 



A^(r) 



he' 



2Vlu 



: Y. %r 



k,i 






>c = clkl 



' ^k.i 



en , = 5r,s ; en ^ • k = ; r, s = 1, 2 . 



{"k,^'4',«} ^^k-k'*^-^ {Ok.r'^k'.s} =0 
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In the above the index k varies from 1 to Nf, k takes 
on all values such that A^(f) is periodic over cube-like 
volumes V and curly brackets imply commutation. 

We now go over to our case of a two dimensional elec- 
tron system by first identifying m to be the effective 
mass m* and introducing the dielectric constant eo into 
the electron-electron interactions. The interactions are 
also modified to take into account the periodic bound- 
ary conditions used (in section II). A is chosen to cor- 
respond to the external magnetic field in the z direc- 
tion i.e. A(f) = {A{f),0). We then introduce approx- 
imations/assumptions into Hq that result in the loss of 
the indistinguishability of electrons. The first of these 
involves the external potential (f>. It is assumed that 
e(j){r) = Vc{z) + Vp{r) where Vp{f) is due to the pinned 
electrons and Vc{z) is assumed to confine the electrons 
at particular values of z. We then assume that Vc{z) 
takes on two different functional forms - the first allow- 
ing for only one state in which the electron is confined 
to the x-y plane and the second allowing for an addi- 
tional state in which the electron is confined a distance 
zo away from the x-y plane (corresponding to the z po- 
sition of the holes bound to the acceptor atoms) . The 
first functional form is assumed to apply to all the elec- 
trons except the recombining electron (index k = c) for 
which the second functional form is assumed to apply. 
This distinguishes the the recombining electron from the 
other electrons. We then make the harmonic approxima- 
tion (as described in section II) which distinguishes all 
the electrons and results in the construction of the two 
Hamiltonians Ti'l and Ti'i corresponding respectively to 
the states in which the z degree of freedom of the recom- 
bining electron is confined to the x-y plane or at a dis- 
tance zq from this plane. As this approximation is made 
only for the x-y degrees of freedom, the z dependence of 
Vp{r) and the electron-electron interactions enters only 
as parameters of the harmonically approximated Hamil- 
tonians Ti-i f. AH the above approximations/assumptions 
are encoded in the following declaration for the approx- 
imated form (Tig'") of the Hamiltonian TYq in which the 
electronic z degrees of freedom are decoupled from the x- 
y degrees of freedom. With the operators written in the 
space of x-y electronic degrees of freedom (g) z electronic 
degrees of freedom, Tip'" may be written as: 



n, = H'iz„p,j 



f 

J2H{zk,PzJ 



where, as mentioned previously, H{z,pz) admits only one 
state localized around z = whereas H'{z,pz) admits an 
additional state localized around z = zg. Pq {Pi) is the 
projection operator onto the multiparticle eigenfunction 
for the z degrees of freedom of the electrons \ijJq) (IV'fo)) 
that has the recombining electron localized around z = 
(z = Zo) - with the corresponding eigenvalue with re- 
spect to Tiz being {—E^). E^ can be considered to be 



the electron-hole binding energy that is an undetermined 
parameter in this theory. 

The approximate form of ?io is now written as: 






(C2) 



which may be seen to have the following eigenfunctions 
that constitute an orthonormal basis: 

{\^i) = 1^,) ® \ro) ® m , i*i=^> = m ® \rzj ® i^^)} 

where (as defined in section II) \ipi) (|?A/)) is an eigenstate 
of 7i^ {Ti.'i) with eigenvalue Ei (Ef) and \iJ;P) is an eigen- 
state of TYq that has a definite number of photons n^ ^ 

in each mode (k, r) (=> eigenvalue E^ = J2]i r ^'^k'^k r)- 
Therefore: 



H?,\^i)^Ei\^i) , Ei- 



E^- 
E' 



EP, 
-EP 



We now approximate Hi (Eq. CI) using: AP{rk) ~ 
y^P(f2) where f^ is the average of the classical equilibrium 
positions of electron k before and after the recombination 
process. This is the dipole approximation_that may be 
justified using input from the experiments that we are 
attempting to model here - the PL peak is at a wave- 
length « 8230A (with the width being very small:« 20A) 
which would be much larger than the spread in the rt 
as computed using |'I'/,_f) or the change in the classical 
equilibrium value of r^ due to the recombination process. 
Therefore: 



771 C ^ — ^ L 

h 



Vk--A{r^) +APM)V 



2m*c^ 



(C3) 



where, in the first summation, the x-y components has 
been separated from the z component of the vectors. 

We now construct the approximate form of 7i (Eq. CI) 
using Eqs. C2 and C3: 



ri — /Tq -|- rij 



(C4) 



7i° is the Hamiltonian in this study. First order time de- 
pendent perturbation theory (Fermi's golden rule) may 
now be applied to compute the intensity of photons emit- 
ted I+{i-o) (or absorbed /-(w)). The computation is now 
analogous to calculations of radiative transitions in atoms 
- here all the Nf electrons constitute the "atom". The 
initial states are declared to be of the form |^/) (recom- 
bining electron localized around z = 0) and the final 
states of the form |\E'f) (recombining electron localized 
around z = zq.) For emission calculations it is assumed 
that that the initial state |^/)+ has no photons and the 
corresponding final state |^f)+ has a single photon with 
a frequency between uj and to + duj whereas for absorp- 
tion calculations the opposite is assumed. Initial states 
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|*/)± that satisfy the requisite constraints on the photon 
numbers are further assumed to be thermaUy distributed. 
The appUcation of the Fermi's golden rule now gives: 



I±{uj)duj 



27r / 



^^'y\±{'^F\n'\\^i)±U{Ei - E, 



F 



Y^^-m 



It may be seen that only one term (— eAf (f^)p^^/?7i*c) 
in Eq. C3, for TiJ, can have a non-zero matrix element 
between the states |^/)± and |^f)± since these differ in 
the z state of the reconibining electron. No other term 
in Eq. C3 has any dependence on the operators corre- 
sponding to the z degree of freedom of the recombining 
electron [zc,'Pz^- ^±(^) is therefore given by the follow- 
ing expression in which the photon states \t\}'^) do not 
appear: 



Af2 






She 



emjPzM)\ 



nm c 



(C5) 



where Afl is the solid angle over which photons of mo- 
mentum k = [uj / c)f are detected (f being a unit vector), 
( )ao is the average over the solid angle, e^ is the z com- 
ponent of e and V is given by Eq. 11 (section II). From 
the above it can be seen that /-(w) = /+(— w). Hence, 
both emission and absorption intensities can be given by 
7+ with the convention that negative w corresponds to 
absorption of photons. An important prediction from the 
above formula is that photons are fully polarized. This 
is most easily seen by noting that the transition involves 
motion of an electron in the z direction, which can couple 
only to the electric fields with a z component. Since the 
polarizations of the photons in the equations below CI 
may be specified with one polarization in the x-y plane, 
this photon does not couple to the transition. Thus, the 
polarization of an emitted photon will be in the k x [^ x k] 
direction. Lastly, using again experimental input that the 
PL peak width is much smaller than the average photon 
frequency tOav, C{uj) may be replaced by C{ujav)- This 
gives: 

I{iu) = C{ujav){{el,)^)AnAnr{uj - h-'E") . (C6) 

This justifies Eqs. 16 and 11 which have E'' — 0. It 
can now be seen that the calculated peak in I(uj) can be 
shifted along the w axis to agree with the experimentally 
observed peak position through a suitable choice of E^. 



APPENDIX D: APPROXIMATIONS AND 
CONVERGENCE CRITERIA 

Here we describe approximations to V{uj) that provides 
a cutoff to the integral in Eq. 12 and also discretizes it 
to a summation. 



It must be noted that P{ll!) is calculated for a finite 
system and therefore it consists of a series of delta func- 
tions as can be seen in Eq. 11. The summations in this 
equation do not go over to integrals as would be the case 
for an infinite system. Therefore, for V(uj) to be a contin- 
uous curve even in the case of a finite system, the delta 
function in Eq. 11 must be replaced by a function with 
non-zero width. We have chosen this function to be a 
Gaussian with width Fq, i.e. 



with g(T,Lu) 



1 



2^ 



27rF 



J2e-^'''Y.\('Pf\^^)\'9i^o,n-\E,^Ef)-cu) 
VaH = ^— 



E' 



-PEi 



The choice of Fq varies with the cases studied - the cor- 
responding values are specified in section III. It can be 
shown that the approximated V{uj) {Va{^)) can be ob- 
tained from V(t) as given by Eq. 13 using: 



Ztt 



+ 00 



dte 



-iujt. 



-Tlt^ 



"2~7'(i) . 



The integral can therefore be cutoff between the limits 
[— T, T] such that the error in Va{i^) due to this {SVa{^)) 
is a small fraction (/) of the expected maximum in Vai^) 
(peak height). Assuming that ■Pa(w) consists of a single 
peak that is Gaussian in shape (this is roughly the experi- 
mental line shapela), the expected peak height in terms of 
the estimated (upper bound) peak width a^^^^ is equal to 
r{t = 0)/\/2^cr^^t. Further it can be shown that \r{t)\ 
is maximum at t = (=1 from Eq. 13). Therefore T is 
chosen such that: 

\6Vaiio)\ < V{0) X 2 Ce'^^dt < ^^4- 

J T y^T^^est 



T|F(1,^)<:^^^ 
' 4' 2 ^- a" 



(Dl) 



(The gamma function F(a, 6) has been defined in Ap- 
pendix A.) / is chosen to be 10~* throughout this study. 
The choice of a"j,i varies with the cases studied - its val- 
ues chosen conservatively to be larger than the expected 
variance of a; in the spectrum being computed. 7'a(w) is 
therefore given by: 



ZTT 



C 

-T 



27r 



with Lo = n—— , n G integers, 

where the discretization of uj follows from the finiteness 
of the time domain integral. 
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We now present a modification to P{t) that aids in 
approximating the above integral to a summation and 
also justifies the assumption that ( the modified) Va{^) 
consists of only one peak. As mentioned in the previous 
section, the parameter A/ in Tif is chosen to be large 
so as to "suppress" the x-y degrees of freedom of the 
recombining electron into becoming "irrelevant" in Tif. 
We have chosen it to be 10'^ e^n^'^/eo throughout this 
study. This makes two normal modes of Ti.'l have very 
large frequencies. These "A modes" may be identified as 
those whose frequency has the strongest dependence on 
the value of A/. It is found (in systems with Nf — 9) that 
the initial states \ipi) do have a significant overlap with 
final states \ipf) whose A modes are excited i.e. then final 
state \ipf) has one or more of the A mode phonons. As 
a result the Va{^) consists of peaks that are separated 
from each other by the A mode frequencies with only a 
small fraction of the total area under the first peak clos- 
est to w = (total area= V^t = 0) = 1). All the peaks 
are centered around points with w < since the ground 
state energy of Tih is greater than that of H.^ (which is 
primarily due to the zero point energy of the A modes in 
Tih as well as because e^ > e™). In other words, in the 
absence of an additional term in Tih corresponding to the 
(negative) binding energy of the electron to the hole all 
the peaks correspond to absorption rather than emission 
of photons (see Appendix C) . We presume that only the 
first peak is relevant to the experiment i.e. we assume 
that the electron-hole binding energy (which is not de- 
termined within our model) is sufficient to shift the peaks 
such that only the first peak comes into the positive lo 
domain. This identification is motivated by the assump- 
tion that the final states \ipf) with no A mode phonons 
best represent the "true" state of the recombined electron 
(with a more realistic Hamiltonian) since, of the different 
\ipf) with varying number of A mode phonons, it is these 
states that maximally localize the recombined electron. 

Having argued that only the first peak is of inter- 
est, the other peaks may by removed from Vaiuj) and 
thereby reduce its width to that of the first peak. This 
is achieved by exponentially suppressing the irrelevant 
peaks i.e. Va{(^) is modified to 7^om(w) given by: 



with 



Vami^) = 



E^ 



-/3-Bi 



E 



e-'^^'x 



E,K^/I^^)I'^"^^''^ + '''^9iTo, h-\E, - Ef) ~ CO) , 

where n^ and rij are the numbers of the two A mode 
phonons in |^/). It can be shown (using arguments simi- 
lar to those leading to Eqs. 12-15) that its Fourier trans- 
form Vamit) is given by: 



-rh^ 



Vamit) =e^r 



oiTr 



Tr 



-pnl 



(D2) 



-l^h 



ih-'n'}t ih-'n';{t + ihp) 



n)^n){^f 



t 



O/-^^), 



a/4 



Vti 



a/ 

„t 



rif = 



^i 



.D, 



(D3) 



where D^ is a diagonal matrix with 7 corresponding to 
the two A modes and zero otherwise. Calculation of the 
traces appearing in Eq. D2 has been outlined in Ap- 
pendix E£3 Hence Vam{^) given by: 



'Pami^L') — ^ 



dte-'^^Vamit) , 



T 



(D4) 



where the cutoff T continues to be given by Eq. Dl. 
The value of 7 is chosen to be 10 throughout this study. 
A larger value could not be chosen since numerical ac- 
curacy (with double precision numbers) of the matrix 
operations needed to compute Vam{t) rapidly decreased 
with increasing 7. At 7 = 10, Vamit) could be cal- 
culated with a relative accuracy of ~ 10~^. In a sys- 
tem with Nf = 9, it has been observed that propaga- 
tion of this error leads to 7-'am('^) having a relative ac- 
curacy of ~ 10"'' even when the limits of integration 

are [-10T, lOT] (T « 61.7 [e^n'Veom*] "^^^ for this sys- 
tem). (Errors are computed relative to a quadruple pre- 
cision calculation.) As this range of integration was the 
largest in this study, we assume that this numerical error 
is negligible in all cases. 

We can now discretize Eq. D4. To do this correctly 
with a minimal number of Vam, it) evaluations i. e. with 
the discretization time being the largest possible, the 
peak of interest must be at the origin. Hence, the quan- 
tity calculated is Vsi^) — Vami^ + LOp), where ujp the is 
expected peak position in Vami^)- In a first approxima- 
tion this is roughly ujq = (-Ef — E^f)/Ti (< 0) where £'; 



'^J 



is the ground state energy of T-C^ r. It is found that tOp is 
actually lower than this because the peak is dominated 
by the overlap of the ground state of Ti^ with excited 
states of Ti 



f 



is therefore calculated assuming that it 
is equal to the expectation value of uj in the distribution 

Vamiuj): 

with ujp ^ ujQ + LU , 

eT-ef Trin,-nf) 
where ujq = —^!- + '— , 

and cZ, = 1 le-*'^o*^™(i) 
I at 



t = 



Vsioj)^^! dte-'i^ + ^P^^Vamit), 

Zn 



(D5) 



-T 

nn 
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The above integration can be transformed into a sum- 
mation without introducing any error if Vs{lo) is zero 
for \uj\ > ujm- Since it is observed that Vs{lo) is roughly 
Gaussian in shape, we assume that this is indeed the case 
with lum — 5 X (3tTi^) where ct^ is the variance of a; in 
the distribution Vs{io) which is given by: 



,-zcjot-5 



dt^ 



I am\J') 



t = o 



Eq. D4 is therefore discretized at a rate > 2 x a;A//27r:l 



■(^) 



with St = — . L = Int 
T' 

L-l 

n— — L 
WK 



(D6) 



to 



T 



L<n< L-l 



and avoid computing 



to now takes on a finite number of values due to the dis- 
cretization of the time domain. We make use of the 

symmetry Vam{-t) = Vamit) 

Vami—t) for negative values of t. It must be noted that 
the large value of ujm in relation to CTi^ is not forced due 
to integration errors (a^j is itself an overestimate of the 
variance due to the peak of interest since there are contri- 
butions from the other exponentially suppressed peaks) 
but due to the presence of a multivalued function (square 
root) in Vamit) whose correct value has been inferred 
only by assuming that its phase does not vary abruptly 
between neighboring t values. It must also be noted that 
the derivatives needed to calculate w and a^^ are com- 
puted numerically. Finally, 'Pam{'-^) is calculated using: 



Vami^) ='Ps{uJ- [Ujp]) , 



(D7) 



where fwr, 



"pi ~ ^p 
different Vs{uj). 
average, 



Wo in cases where there is no averaging over 
However, in cases that demand such an 



is ujp approximated to the nearest integral 
multiple of Si with CTg^( being chosen conservatively so 
that T is the same for all the cases being averaged over. 
This error is expected to be negligible since it is verified 
that 7?i ^ 3(7^ for each of the spectra (7^am(^)) in the 
average. 

It must be noted that the computational scheme guar- 
entees that Vami^) G real numbers since Vam{—t) has 
not been computed independently but has been equated 



to 



' amy^) 



However it must be that 7^am(w) > 0. 
This condition serves as a good test for the validity of 
the approximations. It has been verified that for all our 
results the largest negative value is always a small frac- 
tion (~ 10~^) of the peak value of Vam{^)- 

Convergence criteria: In all the cases studied here the 
only functions needed to compute the potential energy 
terms (and their derivatives) are the gamma functions 
r(a, b) (see Appendix A). As they have been numerically 



evaluated with relative error ^ 10^'^'' we extend the sum- 
mations required to compute the potential (V in Eq. 1 
as computed using Eq. AT) to achieve the same accu- 
racy. Therefore the convergence criterion to determine 
the (classical) equilibrium configuration i.e. the upper 
bound on the root mean squared value of the force per 
electron at the end of the propagation with Hi is cho- 
sen to be lO"^** natural force units (= e^n/eg). Since 
the parameter A/ in Hf is chosen to be 10"^ times the 
natural scale of the same dimension, the convergence cri- 
terion with Tif is weakened, the corresponding bound 
being 10^^^ natural force units. 



APPENDIX E: TRACE COMPUTATIONS WITH 
HARMONIC HAMILTONIANS 

In this section we discuss how the trace operations ap- 
pearing in Eq. D2 are evaluated: 



Z 



G{t) = Tr 



in~^n'}t ih-'^n'iit + ihp) 



with Tif given by Eq. 7, Tih given by Eq. D3, and the 

operators in Tih being related to those in H'l through 
Eqs. 8, 9 and 10. We mention that the exact evaluation 
of the traces, which we do by adapting theoretical tech- 
niques developed in Ref. E^ in the context of coherent 
state properties, is a central feature of our calculations. 
Z is the usual partition function for non-interacting 
bosons. It may be evaluated in the basis of eigenstates 
of H'^ to be: 



Nf 

n 






with Ef = h- 

where w^j is the diagonal element (phonon frequency) of 
the matrix Vti appearing Ti^ . ,— ■ 

Theoretical results of Dodonov and MankoEZi, that use 
the coherent state basis defined by the operators in Tif , 
are employed in evaluating G. This basis, |^), is defined 
by: 

10 - e^l^lO,) , a,j|0,> ^Q- j^ltoNf , (0,|0,> - 1 . 



where |0i) can be identified to be the normalized vacuum 
(ground) state of 7i,f . ^ is a column vector of complex 
numbers of size Nr, and aj ia a row vector of the aj ■ 
operators. It can be shown that the states |^) are over- 
complete and that the trace of an arbitrary operator O 
can be written as: 
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Tr 



O 



d'^fC^-CC 



r^f 



(C|o|0 . 



G niay therefore be evaluated, using additionally the 
cyclic permutation invariance of the trace operation and 
knowing the action of the exponential of the number op- 
erators (a| aij) on the coherent states, as: 



G{t) = Zoe 
U 



ith-^Ef 



^^.-CCi 



T^f 



1/ = 

with E = eH^-/^^)[^^]/2 



(El) 



e *^ '^f* , /i = E*C , z/ = EC 



Therefore, matrix elements of the "propagator" U corre- 
sponding to H'i in the basis |C) must be evaluated. It 
must be noted that this point of view is consistent only 
if the time variable t that appears within H'l (Eq. D3) 
is treated as a constant parameter. This parameter is 
therefore held fixed at fp 7^ and then the matrix el- 
ement is evaluated for all t. Of interest would be the 
value when t = Iq. (It is assumed that the to = value 
is equal to the limiting value as to — > 0.) In other words, 
Tih is viewed as a time independent Hamiltonian. Fur- 
ther, using Eqs. 8, 9 and 10 it can be seen that H'} is 

a second order polynomial in the operators a^. j and a,- . 
The required matrix elements for the most generic time 
dependent Hamiltonian that is of thc-same form has been 
computed by Dodonov and Mankd£3 in a totally differ- 
ent context - we briefly review their solution and then 
specialize to our case. 

The matrix elements of the propagator corresponding 
to the following Hamiltonian arc evaluated in the coher- 
ent state basis: 



n 

h 



[a at] 



,t 



[s] 



K^ 



K) 



where the matrix [K] , the vector [s] and the number r are 
arbitrary complex functions of time (t). (The subscript i 
under the operators has been dropped.) The correspond- 
ing propagator U (and its inverse U^^) may be defined 
through its equation of motion together with an initial 
condition: 



ih—- = nu , U(t = 0)=I 
at 



(E2) 



It can be shown that an invertable operator can be 
defined, up to a multiplicative constant, by its action 
on all the canonically conjugate operators, i.e. if two 
invertable operators U^ and Uy satisfy 



U 



x,{y) 



u, 



x,(y) 



(E3) 



then it must be that lA^ — dAy, where c is a (non zero) 
c-number. This is now used to evaluate the propagator U 
(or equivalently its matrix elements). The multiplicative 
freedom in lA is finally lifted using its equation of motion 
together with the initial condition (Eq. El). Note that 
the elements of b need not be the Hermitian conjugates 
of those in b as the operator U need not be unitary (for 
it could correspond to a non- Hermitian Hamiltonian). 

The above procedure to determine lA can be imple- 
mented only if the operators b and b are known. They 
can be determined from their equation of motion together 
with the initial condition which are obtained using Eqs. 
E2 and E3: 



ifi 



d 


fhl 


in, 


fhl 


1 


fhl 


It 


b 


^[h. 


b 


}• 


b 



(i = 0) = 



(E4) 



where the curly brackets now represent the commutator. 
The above equations are solved by making the following 
anzatz which is motivated by Ti. being a second order 
polynomial in the operators a^ and aj: 



[A] 



[A] 



(E5) 



where A (A) is a c-number matrix (vector) function of 
t. Substituting the above into Eq. E4 it can be seen 
that the anzatz is successful as it leads to the following 
differential equations with appropriate initial conditions 
(the linear independence of the operators a^, a] and 1 
has been used): 



i— = -ASK , Ait = 0) 

dt 



I — = — ASs 
dt 



I 



I 
I 



A(< = 0) = . 



(E6) 



The above may be solved for A and A. An important 
property of the solution is that the matrix A is canonical 
for all t: 



A^SA 



(E7) 



Eq. E3 may now be used to evaluate the matrix el- 
ements of U in the basis |C) after substituting for the 
operators [b] and [b] using Eq. E5 with A and A be- 
ing the solutions of Eq. E6. The matrix elements of Eq. 
E3 (after multiplying from the right by U) in this basis 
reads, with 



Ai A2 
A3 A4 



Ai 
A2 



, {f^\Uit)\u) = U{fi*,u,t) 



|ai|-^ + A2/X* + All U{fj,*,u, t) = uU{ti*,u, t) , 
|a3|-^ + A4/X* + Aaj C/(/x*, i^, t) - ^U{p*,u, i 
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The above can be solved (assuming the existence of Aj~ 
and using Eq. E7) for the /x* and u dependence of 
U {/J,* , u , t) . The unsolved time dependence appears as 
an arbitrary multiphcative integration constant: 

\n[U{ti*,u,t)] = ^[-,1^ A^'A2fi* + ly'^AsA^'u] 
+ti^A^^ [v - Ai] + v^ [\2 - AsAr'Ai] + /(t) . (E8) 

Using the above and evaluating the matrix elements of 
Eq. E2 (for U), it may be seen (using Eqs. E6 and E7) 
that Eq. E2 is indeed satisfied in its /x* and v depen- 
dences, and it leads to the following equations determin- 
ing /: 



dt 
1, 



-Tr 



Id 
2dt 

dA, 
^ dt 



AfAgA^^Ai 



(E9) 



dX^ 
dt 



-Ai 



fit - 0) - 



Hence, as argued previously, Eq. E3 leaves only a c- 
number multiplicative freedom in U which is lifted by 
using Eq. E2. This completes the generic solution of the 
matrix elements of interest. 

We now specialize to our case to evaluate G (Eq. E4) 
by identifying H to be H'i. This implies (using Eqs. 8, 
9, 10 and D3): 



G{t) (Eq. El) is now determined as the matrix ele- 
ment (/x|Z^|i/) = U{fj.*,i',t) (usiB§ Eq. ElO) is given by 
Eq. E8. (Dodonov and MankocZl also obtain this ma- 
trix element for the special case of ftf ex I and D-, = 0) 
Further, the integral in Eq. El may be carried out an- 
alytically (assuming the positive definiteness of the ap- 
propriate matrix) to give: 



In 



Git) 



Zn 



(Ell) 



ih-'Ef + fit) 

ln{det[AS+]} lj,T,-ip 
2 2 



A = S I EJE , F = El 



sre S+ = 


I 

1 


, A = S 


with J = 


" A-IA2 A-i 
[A-^f AsAr^ 


and 1 — 


A2 


Ai 

AsAr'Ai . 



E 



E 
E 



This completes the analytical determination of the re- 
quired trace operations. 






K^T'flfT, s^T'flfW, r = 



h-h"" 



f ■ 



The above may now be used to obtain explicit solutions 
to Eqs. E6 and E9 - as mentioned previously the time 
t that appears within Ctj is treated as a constant and 
then the solution is evaluated for t being equal to this 
constant. This gives, with 

p^jifi/t + D^] ^ p^j-inft~-D^] ^ 

and knowing (using Eqs. 9 and Bi) that T is canonical, 
i.e. T^ST = S: 



Ai(t) = UtpU- V'^PV* , 
A2(t) = UtpV- V^PU* , 
Asit) = U^PV*- Vtpu , 
Aiit) = U^PU* - Vtpv , 

Ai(t) = Ut [P-I] A- V^ 



(ElO) 



Mt) = u^ 



Vt [P - II A 



1 ^T, 



fit) = - [A^AaAr'Ai - In (detAi) 



+ — Aj^Ai — A2A2 2A2A]^ 

- A^Sinh [D^] {2Cos [f7/i] - Cosh [D.^]} A 
-iA^^Cosh [D.^] Sin [^ft] A - ih~\ft . 
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